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Magnitude-yield  calibration  of  Northern  Eurasia  can  be  improved  by  using  the  large  chemical  and  Peaceful  Nuclear 
Explosion  (PNEs)  seismic  datasets  acquired  by  the  Soviet  Deep  Seismic  Sounding  (DSS)  program.  Currently,  we 
are  using  inversion  and  numerical  modeling  to  analyze  the  coda  amplitude  parameters  and  P-wavc  travel  times  from 
numerous  DSS  explosions,  which  would  lead  to  a  detailed  empirical  mapping  of  seismic  properties  within  northern 
Eurasia. 

In  this  report,  PNE  Lg  coda  amplitudes  from  seven  profdes  are  inverted  to  produce  maps  of  the  following  new  coda 
parameters:  1)  frequency-dependent  scattering  intensities,  2)  geometrical  parameter  (y),  and  3)  effective  quality  (Qe). 
The  rand  frequency-independent  Qe  parameters  replace  the  traditional  frequency-dependent  Q(f)  =  Qoiflfo)'1 
parameterization  and  provide  stability,  independence  of  assumptions  of  a  background  model,  and  convenience  in 
interpretation.  As  we  have  shown  before,  parameter  y  also  yields  itself  to  independent  numerical  waveform 
modeling,  correlates  with  crustal  tectonic  types,  and  may  be  a  portable  attribute  valuable  for  coda  regionalization. 

All  three  maps  show  notable  correlations  with  regional  tectonics  and  with  recently  derived  maps  of  phase  amplitude 
ratios  and  Lg  Q. 
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OBJECTIVES 

The  correct  and  transportable  measurement  of  attenuation  parameters  is  critical  in  coda  magnitude  studies.  Recently 
(Morozov  et  al.,  in  press)  we  showed  that  short -period  coda  attenuation  of  regional  arrivals  can  be  interpreted  in 
terms  of  the  geometrical  spreading  (y)  and  effective  attenuation  (Qe)  parameters.  Taken  together,  they  represent  the 
generalized  attenuation  coefficient:  u(J)  =  y+  7rflQe.  This  model  removes  the  commonly  observed  ambiguity  of  the 
frequency-dependent  Q{f),  and  parameter  /becomes  most  useful  for  interpretation.  Modeling  of  P- wave  coda  at 
regional  distances  and  analysis  of  two  Peaceful  Nuclear  Explosions  in  Russia  in  terms  of  this  concept  showed  a 
remarkable  agreement,  particularly  in  the  values  of  y . 

In  this  report,  we  apply  the  same  parameterization  to  the  inversion  of  spatially-variable  surface-consistent  Lg  coda 
parameter  cdj)  from  five  PNE  profiles  in  Siberia.  Minimalistic  model  parameterization  and  the  Constrained 
Simultaneous  Iterative  Reconstruction  Technique  (SIRT)  are  used  for  unbiased  inversion.  The  resulting  map  of  the 
geometrical  parameter  (y)  is  in  excellent  agreement  with  both  modeling  and  inversion  in  Morozov  et  al.  (in  press). 
The  values  of  /differentiate  between  the  crustal  tectonic  types  within  the  study  area  and  are  also  in  agreement  with  a 
summary  of  over  forty  worldwide  studies  (Morozov,  in  press).  In  addition,  scattering  amplitudes  are  found  to  be 
highly  variable  and  strongly  correlate  with  geological  structures.  Coda  Qe~l  values  also  show  pronounced 
correlations  with  geological  structures,  although  these  values  may  also  be  somewhat  overestimated  compared  to  two 
detailed  PNE  point  studies  by  Morozov  et  al.  (in  press). 

RESEARCH  ACCOMPLISHED 


Rationale 

Geophysical  characterization  of  the  Earth’s  lithosphere  broadly  consists  of  descriptions  of  its  velocity,  density,  and 
attenuation  structures.  Combinations  of  these  parameters  are  used  in  a  broad  range  of  applications,  from  providing 
structural  and  rheological  constraints  for  geodynamics  to  event  location,  discrimination  of  seismic  events, 
calibration  and  regionalization.  For  the  first  two  of  these  parameters,  a  variety  of  precise  methods  have  been 
developed,  but  in  situ  attenuation  remains  notoriously  difficult  to  measure  and  interpret.  In  particular,  while  quality 
parameter  ( Q )  often  showing  strong  frequency  dependence,  characterization  of  such  frequency-dependent 
attenuation  (Q(f))  may  be  ambiguous  and  dependent  on  model  assumptions  and  inversion  methods. 

Morozov  (in  press)  recently  proposed  a  new  technique  for  measuring  attenuation  without  the  use  of  the  ambiguous 
Q(j)  and  illustrated  it  on  a  number  of  examples.  In  this  study,  we  employ  this  method  to  invert  coda  Q  values  from 
five  Peaceful  Nuclear  Explosion  (PNE)  profiles  in  Russia  and  to  map  the  frequency-dependent  attenuation  across  a 
broad  area  of  Siberia  (Figure  1). 

Different  types  of  attenuation  (body,  surface,  Lg,  and  coda  waves)  are  often  presented  in  the  following  power-law 
form  of  frequency-dependence: 

e(/)=a(/0’  (1) 

where  Q  is  the  apparent  quality  parameter,  and  /fl  is  some  reference  frequency  often  taken  to  equal  1  Hz  (Aki,  1980). 
However,  as  argued  by  Morozov  (in  press)  and  Morozov  et  al.  (in  press),  parameters  of  such  a  power-law 
dependence  may  be  strongly  influenced  by  the  observation  process,  and  in  particular  by  the  assumptions  made  about 
the  geometrical  spreading  of  scattered  (in  the  case  of  coda  Q )  or  other  waves.  Parameters  Qo  and  ij  trade  off  with 
each  other  through  the  uncertainty  of  the  assumed  geometrical  spreading.  The  geometrical  spreading  (GS),  in  its 
turn,  cannot  be  considered  constant  within  large  areas,  and  -10%  uncertainties  in  it  can  account  for  much  of  the 
observed  frequency  dependence  of  the  apparent  Q  (Morozov,  in  press).  Unconstrained  variations  of  GS  also  make  it 
difficult  to  correlate  attenuation  parameters  from  different  regions,  and  also  to  link  them  to  lithospheric  properties. 

In  seismic  calibration  and  nuclear  test  monitoring  studies,  uncertainties  in  Q0,  rj.  and  geometrical  spreading  are 
particularly  undesirable,  because  they  prevent  the  construction  of  consistent,  seamless  models  of  large  areas,  and 
complicate  magnitude  calibration. 
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Numerical  waveform  coda  modeling  (Morozov 
et  al.,  in  press)  showed  that  for  a  layered  1-D 
model  with  a  constant,  frequency-independent 
crustal  attenuation,  the  observed  coda  Q{f) 
exhibited  a  strong  dependence  on  frequency. 

Although  such  dependence  is  often  observed  in 
coda  Q  and  Lg  Q  data  (e.g.,  Der  et  ah,  1986; 

Campillo,  1987,  1990;  Frankel,  1990;  Mitchell 
and  Cong,  1998;  Benz  et  ah,  1997;  Mitchell  et 
ah,  1997;  McNamara  et  ah  1996;  McNamara, 

2000;  and  Erickson  et  ah,  2005),  this  modeling 
shows  that  it  may  relate  to  the  crustal  structure 
(particularly  to  layering)  and  not  necessarily  to 
its  rheology,  fracturing,  or  fluid  content.  Similar 
observations  were  made  earlier  by  Anderson  et 
ah  (1965)  for  long-period  surface  waves  and  by 
Mitchell  (1991)  for  Lg  waves  in  the  Basin  and 
Range.  Strictly  speaking,  because  the  values  of 
Q  and  77  trade  off  with  the  imprecisely  known 
velocity/density  structures,  these  parameters  are 
not  necessarily  be  associated  with  attenuation 
properties  and  may  not  reliably  correlate  with 
geological  structures  or  with  the  variations  of 
the  physical  state  of  the  lithosphere.  However, 
the  potential  association  of  the  widespread 
observations  of  frequency-dependent  Q(f)  with  variations  of  the  structure  still  has  not  been  examined  ,  and  values  of 
Q  are  typically  interpreted  as  related  to  rheology,  physical  state,  or  fluid  content  of  the  lithosphere  (e.g.,  Mitchell  et 
al.,  1997). 

In  the  practice  of  Q(f)  interpretations,  and  in  particular  in  Lg  coda  Q  studies  (such  as  coda  magnitude  calibration), 
the  sensitivity  of  coda  Q  to  the  frequency  is  usually  overcome  by  referencing  the  coda  Q(J)  dependence  to  some 
common  frequency  level,  such  as  1  Hz,  in  all  comparisons  (e.g.,  Mitchell  et  al.  1997;  Phillips  et  al.,  2004).  However, 
this  still  does  not  resolve  the  problem  of  Q(f)  ambiguity,  as  such  referencing  process  depends  on  the  velocity 
structure  through  parameter  77,  making  the  resulting  values  less  portable.  Note  that  because  of  the  inherent 
dependence  of  the  results  on  unconstrained  GS,  changes  in  inversion  approaches  sometimes  also  lead  to  updated 
values  of  Q0  and  77  (for  example,  indicated  by  Xie  et  al.  (2006)  in  relation  to  Mitchell  et  al.  (1997)  or  Xie  (2002)  in 
relation  to  McNamara  et  al.,  (1996)).  Thus,  the  interplay  of  the  assumed  forms  of  GS  and  power-law  Q(f)  in 
dependence  (1)  creates  a  difficult  environment  for  measurements  and  interpretation  of  attenuation. 

in  many  coda  Q  and  Lg  Q(f),  and  particularly  in  nuclear  test  monitoring  studies,  Q  is  treated  as  a  sort  of  coda-shape 
or  spectral-amplitude  shape  parameter.  When  used  strictly  for  matching  the  observed  time-frequency  dependent 
amplitudes,  the  above  model  works  even  with  imprecise  or  even  arbitrary  GS.  However,  the  ability  of  matching  the 
amplitudes  still  does  not  prove  that  the  propagation  model  is  correct  and  can  be  used  to  constrain  the  crustal  and 
mantle  properties.  This  problem  is  well-known  in  the  inversion  theory  -  under-constrained  (or  over-parameterized) 
models  tend  to  fit  the  data  perfectly  while  leaving  ambiguities  in  the  solution  (Menke,  1989).  When  interpreting 
such  solutions  and  comparing  them  to  each  other,  it  is  important  to  avoid  using  parameter  combinations  that  are 
related  to  the  unconstrained  degrees  of  freedom.  Unfortunately,  both  Q0  and  77  in  eq.  (1)  are  such  parameters,  which 
co-vary  with  the  unconstrained  GS. 

To  resolve  this  difficulty  of  intermixed  GS  and  attenuation  effects  and  the  resulting  apparent  character  of  Q(f,  we 
proposed  (Morozov,  in  press,  and  Morozov  et  al.,  in  press)  to  incorporate  a  geometrical  spreading  parameter  y  in  the 
following  equation  for  coda  amplitude  decay: 

4,/)=  AXfY'e-^  =  (2) 

where  k=  n!Qe,  and  vis  the  theoretical  GS  exponent  corresponding  to  the  type  of  scattered  wave  modes  that  is 
considered  as  predominant.  Arguments  for  form  (2)  based  on  the  scattering  theory  were  given  by  Morozov 


Figure  1.  Sources  (red  circles)  and  receivers  (blue  triangles) 
with  measured  frequency-dependent  coda 
</ln(Amplitude)/df  derivatives.  Black  stars  are  the 
PNEs,  and  small  black  triangles  -  3-component 
receiver  locations.  Only  3°  grid  cells  containing 
sources  or  receivers  were  used  in  the  inversion.  Two 
PNEs  used  in  Part  I  are  indicated:  Q4  -  Quartz-4, 
and  K3  -  Kimberlite-3  (see  also  Figure  2). 
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(unpublished);  however,  this  expression  can  also  be  viewed  as  a  purely  empirical  approximation  (Morozov  et  al., 
2006;  Morozov,  in  press).  The  subscript  “e”  in  Qe  indicates  the  “effective”  character  of  the  quality  parameter,  in 
which  we  do  not  differentiate  between  the  intrinsic  attenuation  and  small-scale  (Rayleigh-  and  Mie-type)  scattering. 
Note  that  parameter  k ,  and  consequently  also  Qe,  can  be  frequency-dependent.  However,  no  indications  of  such 
dependence  were  found  in  several  published  examples  of  Lg,  and  body -wave  data  (Morozov,  in  press),  and  also  in 
surface-wave  and  normal-oscillation  data  (Morozov,  unpublished). 

In  parameterization  (2),  the  two  parameters  ( Q0 ,  rj)  of  the  power  law  (1)  are  replaced  with  three:  v,  y and  k.  Within 
the  typical  errors  of  spectral  observations,  there  is  little  room  for  seeking  frequency-dependent  y  or  k,  and  they  can 
therefore  be  considered  frequency-independent.  Morozov  (in  press  and  unpublished)  re-examined  several  surface-, 
body-,  Lg-,  and  coda-wave  studies  and  showed  that  model  (2)  was  invariably  applicable  to  all  of  them.  In  most 
cases,  this  resulted  in  significant  increases  of  the  reported  values  of  Q. 

Parameter  v  significantly  trades  off  with  y,  and  therefore  it  is  acceptable  to  use  amplitudes  A(t,f)  corrected  for 
geometrical  spreading  with  fixed  v values  determined  from  theoretical  considerations,  as  it  is  commonly  done  in  Q(f) 
studies.  In  such  cases,  parameter  /describes  the  spatial  variations  of  uncompensated  geometrical  spreading.  In  the 
case  of  regional  Lg  coda  Q  of  this  study,  we  use  a  model  of  cylindrically-propagating  coda  waves  generated  within 
the  near-surface.  In  this  case,  the  spreading  area  of  scattering  surrounding  the  source  and/or  receiver  compensates 
the  geometrical  spreading  of  coda  waves,  and  vcan  be  set  equal  to  0  (Morozov  and  Smithson,  2000).  However,  this 
is  only  a  simplified  model,  and,  for  example,  diving  and  reflected  scattered  waves,  deeper,  or  multiple  scattering 
would  change  the  f  1  law,  and  make  it  regionally  variable.  Our  objective  in  this  study  is  to  measure  these 
variations,  in  the  form  of  a  spatially-variant  y. 

Parameter  /was  found  to  correlate  with  lithospheric  types  and  ages;  in  particular,  the  value  of  yD  =  0.008  s"1 
appeared  to  clearly  separate  the  tectonically  stable  (/<  yn)  and  active  areas  (/  >  /n;  Morozov,  in  press).  In  addition, 
numerical  modeling  of  the  regional  P-wave  coda  conducted  in  several  PNE-based  models  (Morozov  et  al.,  in  press) 
showed  that  the  modeled  /was  in  a  close  agreement  with  the  measurements  from  PNE  data,  and  coda  Q  could  be 
related  to  the  S-wave  crustal  Qs.  Therefore,  it  is  particularly  interesting  to  look  for  correlations  of  /  with  geological 
structures  across  a  large  contiguous  area  like  the  PNE  profiling  area  of  this  study. 

From  the  viewpoint  of  the  classical  scattering  theory  (e.g.,  Chernov,  1960,  p.53;  Nikolayev,  1968),  parameterization 
(2)  is  merely  a  return  to  the  traditional  frequency-dependent  attenuation  coefficient  a{f)  instead  of  the  “scattering 
quality”  parameter  Qs  introduced  by  Aki  (1980).  Dainty  (1981)  gave  a  justification  for  expression  (2)  as  describing 
scattering  of  1-30-Hz  5-wave  on  larger  scatterers.  In  coda  Q  studies,  such  scattering  could  occur  by  means  of  diving, 
reflected,  and  mode-converted  waves  propagating  through  the  crust.  Note  that  such  processes  mostly  involve 
deterministic  effects,  such  as  ray  bending  in  crustal  velocity  gradients  and  reflections  from  major  lithospheric 
structural  contrasts,  and  therefore  the  use  of  Q,  leads  to  undesirable  connotation  of  random  scattering  and  loss  of 
correlation  with  the  structure.  Therefore,  instead  of  using  Qs,  we  prefer  using  the  uncompensated  GS,  which  simply 
represents  the  frequency-independent  part  of  the  attenuation  coefficient:  /=  a(f=0). 

Preliminary  analysis  of  two  PNEs  (Figure  1;  Morozov  et  al.,  in  press)  suggested  that  /values  could  be  remarkably 
stable,  even  in  the  case  of  surprisingly  strong  frequency  dependence  (  rj «  1)  within  the  Siberian  Craton.  The 
observed  values  of  /also  agreed  well  with  those  derived  from  waveform  synthetics  computed  for  the  study  area.  We 
argued  that  /could  potentially  provide  a  stable  and  transportable  criterion  for  correlating  the  observations  of 
attenuation  with  crustal  tectonic  types. 

In  this  paper,  we  further  apply  the  approach  arising  from  parameterization  (2)  to  mapping  the  short-period 
attenuation  structure  in  central  Siberia  by  using  several  PNE  datasets  (Figure  1).  Detailed  descriptions  of  the  PNE 
profiles  and  geological  setting  of  the  study  area  were  recently  given  by  Li  et  al.  (2007  and  submitted  to  BSSA)  and 
are  not  repeated  here.  The  resulting  / and  Qe  maps  show  good  stability  and  correlations  with  velocity  structures  and 
tectonics.  Because  the  parameters  being  mapped  are  directly  related  to  crustal  and  upper  mantle  properties,  they 
should  also  provide  a  good  basis  for  comparisons  to  other  studies  and  for  incorporation  into  forward  modeling. 

Data 

The  dataset  consists  of  log-amplitude  coda  amplitude  readings  from  the  PNE  explosion  records  from  five  Deep 
Seismic  Sounding  profiles  in  Russia  (Figure  1).  Figure  2  shows  an  example  of  transverse-component  (relative  to  the 
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source-receiver  direction)  record 
from  PNE  Kimberlite-3  (labeled 
K3  in  Figure  1)  located  near  the 
edge  of  the  Siberian  Craton.  Note 
the  high  density  of  recordings  (10- 
15  km  spacing)  and  the 
differences  in  the  wavefield 
propagating  within  and  beneath 
the  West  Siberian  Basin  (west  of 
the  PNE)  and  the  Siberian  Craton 
(Figure  2). This  PNE  was  also 
used  in  detailed  studies  by 
Morozov  et  al.  (in  press). 

The  data  were  carefully  edited  by 
removing  poor  and  clipped 
records,  and  all  regional  arrivals 
were  identified  and  picked.  The 
Lg  distance  ranges  and  time 
windows  were  selected 
interactively  from  all  PNE  record 
sections.  Coda  windows  started  20 
sec  after  the  picked  Lg  onsets  and 
extended  for  50-  100  sec.  Pr e-Lg 
noise  windows  (i.e.,  the  P-  and  S- 
wave  coda;  cf.  Morozov  and  Smithson,  2000)  were  also  picked.  By  using  these  windows,  Lg  to  pre-Zg  amplitude 
ratios  were  calculated,  and  traces  with  these  ratios  below  1.1  were  discarded. 

Edited  records  were  further  band-pass  filtered  within  four  overlapping  frequency  bands  of  1-2,  1.5-3,  2-4,  and 
3-5  Flz.  Within  each  band,  three-component  trace  envelopes  were  formed,  and  d\n(Amplitude)/dt  derivatives  and 
their  standard  error  estimates  were  measured  by  using  the  “robust  fit'’  technique  in  Matlab.  The  resulting  values  of 
frequency-dependent  log-amplitude  coda  slopes  were  saved  in  a  database  which  was  used  in  the  subsequent 
analysis. 

Surface-Consistent  (A,/,  Qe)  Inversion  Method 

The  power  and  utility  of  the  cdj)  model  is  in  its  direct  link  to  the  observations  (Morozov,  in  press).  In  the  case  of 
regional  coda  attenuation,  coda  ac(f)  simply  becomes  the  negative  time  derivative  of  the  logarithm  of  coda 
amplitude  (2): 

(3) 

dt 

which  can  be  measured  by  fitting  the  coda  envelopes  in  each  individual  record.  Because  of  the  source-receiver 
reciprocity,  the  observed  coda  amplitude  (2),  and  consequently  ac(J),  contain  contributions  from  scattering  near  the 
source  and  receiver.  Spatial  separation  of  these  two  areas  and  the  corresponding  differences  in  their  cdj)  values 
could  be  significant  for  mapping,  particularly  in  the  case  of  the  long-range  PNE  sources.  Thus,  we  need  to  invert  the 
observed  ac(f)  for  spatial  variations  in  a(f).  To  accomplish  this,  we  use  the  surface -consistent  approach,  and 
represented  cdj)  as  functions  of  the  coordinates  and  frequency.  The  surface-consistent  model  implies  that  for  the 
collocated  source  and  receiver,  the  values  of  cdj)  should  be  the  same,  and  for  spatially  close  locations,  odj)  values 
are  close. 

To  derive  a  model  for  the  observed  y  and  Qe,  we  use  an  approximation  for  the  event-energy  envelope  by  Morozov 
and  Smithson  (2000),  in  which  the  intensity  (energy  density)  of  the  primary  event  was  described  by  a  short  pulse  of 
parabolic  shape  arriving  at  time  to,  with  amplitude  Pt,  and  duration  r.  The  coda  of  this  arrival  was  described  by 
exponentially  decaying  amplitude  following  the  pulse: 


West-Siberian  Basin  Siberian  Craton 


Figure  2.  Transverse-component  record  from  PNE  Kimberlite-3  (K3  in 
Figure  1),  filtered  within  0.5  -  8.0  Hz  pass-band.  The  regional 
phases  (labeled  P,  Pg,  S,  and  Lg)  are  clear  and  observed  to  far 
source-receiver  distances.  Note  the  difference  between  the 
branches  of  Lg  in  two  directions  from  the  PNE. 
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where  A{f)  is  the  scattering  efficiency  factor  (the  relative  coda  intensity  at  the  time  of  the  primary  event).  In  general, 
this  parameter  may  also  depend  on  the  frequency. 


From  PNE  Quartz-4  records  (Q4  in  Figure  1),  the  duration  of  the  primary  onsets  was  estimated  as  z=  1.25  s  for  the 
teleseismic  P,  and  A  =  0.22  for  all  events  (Morozov  and  Smithson,  2000).  The  relative  coda  intensity  at  the  time  of 
the  primary  event  was  0.27  (Equation  [4]),  which  can  be  considered  as  significantly  smaller  than  1.  Therefore, 
the  total  wavefield  intensity  at  the  receiver: 


PSR  (0  ^  [Pprimary  (0  +  Pcoda  (A  4  »  <*S  )J*  [Ppn,nan<  (0  +  Pcoda  (A  4 


(5) 


(quantities  with  subscripts  S  and  R  correspond  to  the  source  and  receiver  locations,  respectively,  and  the  asterisk 
denotes  time  convolution)  can  be  approximated  to  the  first  order  in  A: 

PSR  (0  —  Pprimary  (0  *  Pprimary  (0  Pprimary  (0  *  \Pcoda  (A  4  ’  as  )  Pcoda  (A  4  >  aR  )]  '  (^) 


The  resulting  coda  intensity  is  given  by  the  second  term  in  this  expression: 

PsRjvdSt) «  Axs ,  » f)e~2a(xs’ys’f)t  +  A(xr , yR  ?  > ,  (7) 


where  we  also  assumed  surface  consistency  of  A.  From  eq.  (2),  u(x,y.J)  is  further  approximated  by  a  linear  function 
of  frequency: 

a{x,  y,  /)  =  y(x,  y)  +  k{x,  y)f ,  (8) 


Taking  into  account  the  source  and  receiver  site  effects,  the  coda  intensity  (7)  is  further  approximated  as: 


'  SR, coda 


-2  act 


(9) 


where  C(f)  is  a  time-independent  amplitude  normalization  factor  of  the  specific  seismic  record.  With  these 
approximations,  eq.  (9)  leads  to  the  following  system  of  non-linear  equations: 


d . .  =  +  -C(/)=0 


(10) 


for  all  t,f  and  record  numbers  n.  Solving  equations  (10)  is  equivalent  to  minimizing  the  following  objective 
function: 


®(r,^,4>= jlIJ[4(/> 


\e  +AR{f)e 


-2[  a(xR,yR,f)-ac\ 


Clf)]dt , 


f  ”  T„ 


(11) 


where  T„  is  the  coda  observation  time  window.  The  summation  in  eq.  (11)  takes  place  over  all  records  and 
frequencies  of  interest,  and  the  minimization  is  performed  in  terms  of  parameters  y,  /c  A(f,  and  C(f).  In  this  study,  to 
improve  the  stability  of  the  inverse,  we  did  not  consider  the  dependence  of  k.  A,  and  g on  frequency,  and  solved  for 
variations  of  y,  k  and  A  as  functions  of  spatial  coordinates  only. 

To  define  the  spatial  inversion  grid  for  parameters  y  k.  and  A,  care  needs  to  be  exercised  in  order  to  not  over¬ 
parameterize  the  model  excessively.  To  achieve  a  “minimalistic”  parameterization,  we  defined  a  spatial  grid  with 
increments  of  3°  in  both  latitudes  and  longitudes  covering  the  entire  area  of  profiling.  Only  cells  containing  the 
actual  source  or  receiver  points  (Figure  1)  were  included  in  the  inversion,  which  resulted  in  53  cells.  For  each  of 
these  selected  grid  cells,  a  bilinear  basis  function  (finite  element)  y/,{x,y)  was  defined  ( x  was  the  longitude,  y  is  the 
latitude,  and  i  is  the  cell  number),  such  that  y/,(x,y)  =  1  in  its  center  and  y/,(x,y)  ~  0  at  the  centers  of  all  other  cells. 
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“Blocky”  basis  functions  typical  in  travel-time  tomography  (with  t//,(xjj  =  1  within  the  entire  cell)  were  also  tried, 
with  no  significant  differences  in  the  results. 

With  properly  defined  basis  functions,  the  following  identity  holds  at  any  point  (x,y): 


Z  V'i(x,y)=  !’ 


(12) 


where  the  summation  includes  all  cells.  With  the  help  of  this  functional  basis,  values  of  y,  k,  and  A  at  any  point  ( x ,  y) 
can  be  expressed  through  the  corresponding  values  (y,  /r„  A,j  at  the  nodes;  for  example,  for  y. 

r{x,y)=J]f//i(x,y)rr  (13) 


Within  the  class  of  functions  given  by  the  finite-element  functional  basis  (13),  the  objective  function  (11)  becomes  a 
non-linear  function  of  model  parameters:  ®(y,  a),  A„  C„).  However,  unconstrained  degrees  of  freedom  are  still 
present  in  this  parameter  space,  such  as  an  arbitrary  scaling  of  parameters  A  and  £  in  eq.  (11)  and  source-receiver 
trade-off  caused  by  non-uniform  sampling  of  the  model.  These  degrees  of  freedom  are  typically  removed  by 
regularization,  which  is  performed  by  adding  constraint  terms  to  the  objective  function  (11),  for  example: 


®(yi,Ki,Ai,£k)=<b(yi,Ki,Ai,£t )+  £ 


(14) 


where  functions  ¥1  (y,  k„  A, )  penalize  various  undesirable  types  of  instabilities.  However,  this  method  of 
regularization  requires  a  non-trivial  selection  of  weights  wy  and  biases  the  solution  from  minimizing  the  objective 
function  (11)  and  satisfying  equations  (10).  Most  importantly,  properties  of  the  resulting  constraints  may  be  difficult 
to  assess  in  order  to  evaluate  their  effects  on  the  interpretation.  Therefore,  a  simpler  regularization  scheme  with 
readily  interpretable  constraints  would  be  preferable. 

In  order  to  obtain  a  stable  solution  without  being  unbiased  from  the  minimum  of  <D(y,  Kh  Ah  <Q,  we  utilized  the 
iterative  nature  of  the  SIRT  (Simultaneous  Iterative  Reconstruction  Technique)  solver  that  we  employed.  After  each 
iteration,  the  solution  was  corrected  to  satisfy  the  following  criteria: 

1)  A,  >  0  and  x;  >  0  at  all  points.  These  constraints  correspond  to  non-negative  scattering  amplitude  and  Q, 
respectively.  Note  that  unlike  xrh  y  can  take  both  negative  and  positive  values. 

2)  ^  Aj  =  NgM  ,  where  Ngrii  is  the  total  number  of  grid  points  involved  in  the  inversion.  This  constraint  was 

i 

enforced  by  rescaling  all  values  of  A,  and  removing  the  scaling  invariance  of  the  minimum  of  expression  (11). 

3)  |/L;  -  /^mooth|  <  /?/i™00th  ,  where  ^mooth  is  the  smoothed  value  of  T,  derived  by  averaging  the  adjacent  points 

(and  excluding  the  i-th  one),  and  [)  is  the  tolerance  parameter  selected  equal  to  0.05  in  our  inversion.  This 
criterion  guaranteed  that  adjacent  cells  differ  by  no  more  than  f:i  in  terms  of  A,-.  Similar  constraints  were 
implemented  for  y  and  zq. 

4)  Additionally,  similar  constraints,  for  example:  | A.  -  Afnon (x,  y)  <  /?'  /}?pnon (x,  y) ,  could  be  added  to  keep  the 

solution  in  the  vicinity  of  some  a  priori  model.  Such  constraints  could  be  useful  to  remove  instabilities  at 
isolated  points  near  the  circumference  of  the  model,  or  for  extrapolation  of  the  solution  outside  of  the  area  of 
coverage.  The  a  priori  model  can  be  parameterized  at  user-specified  sets  of  support  points  between  which  the 
values  of  Aap"on(x,  v)  are  linearly  interpolated  by  using  a  Delaunay  triangulation.  However,  we  did  not  use 
such  constraints  in  the  solution  presented  below. 


To  invert  for  y,  x),  Ah  and  C„  numerically,  we  used  the  SIRT  method  by  applying  it  directly  to  Equation  10.  For  a 
small  parameter  perturbation,  these  equations  become: 

df.n,t  =  d°  +SdfM  =  d°  +  [SAs  -  2 t(Sys  +  f  5ks  +  [SAR  -  2 t(SyR  +  /8k R  -  SC„  =  0.  (15) 
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Quantities  8y,  8k ;  and  8A  evaluated  at  the  source  and  receiver  locations  are  further  expressed  as  combinations  of  the 
corresponding  values  at  the  nearest  grid  nodes  by  using  equation  13.  By  combining  all  parameter  perturbations  Sy, 
Sk„  8Ah  and  8Cn  into  one  model  vector  p,  the  resulting  equations  form  a  linear  system: 

A  dip  -  B  =  0  ,  (16) 

which  was  over-  or  mixed-determined  (Menke,  1989;  7,986  data  points  and  885  unknowns  in  our  inversion).  Its 
approximate  solution  in  the  least-squares  sense  is: 

4)  =  (ArA)_1ArB,  (17) 

where  T denotes  the  matrix  transpose.  In  SIRT  (and  related  back-projection  methods),  the  diagonally-dominant  A rA 
matrix  is  replaced  in  this  equation  with  its  diagonal,  whose  inverse  can  be  easily  calculated: 

dp  *  c [diag (a r A )]~ '  A 7  B  .  (18) 

This  method  requires  the  storage  of  only  two  vectors:  A  B  and  the  diagonal  of  AZA,  and  consequently  it  can  be 
efficiently  applied  to  very  large  problems.  In  our  study,  the  inverse  (18)  was  calculated  during  iterative  scanning 
through  the  entire  sc(f)  dataset,  and  was  terminated  when  model  updates  became  sufficiently  small.  Between  the 
iterations,  “trimming”  operations  were  applied  to  the  model  in  order  to  make  it  satisfy  the  constraints  above.  Our 
final  model  (Figure  3)  resulted  in  ~90%  data  error  reduction  from  the  starting  model  with  y  =  0,  /q  =  0,  and  7,,  =  1 . 

To  suppress  parameter  oscillations  during  this  process,  a  scaling  factor  c  <  1  (c  =  0.5  in  our  inversion)  was  applied 
during  stepping  (18)  and  was  gradually  reduced  when  data  error  increases  were  detected. 

Results 

The  resulting  maps  of  A,  y ,  and  k  within  the  area  of  Lg  coda  Q  data  coverage  from  five  PNE  profiles  are  shown  in 
Figure  3. To  produce  continuous  maps,  we  used  linear  interpolation  within  a  Delaunay  triangulation  constructed  on 
the  centers  of  the  53  grid  cells  that  were  used  in  the  inversion  (Figure  1).  This  explains  the  shapes  of  some  of  the 
features  located  near  the  edges  of  the  coverage  area  (such  as  the  Baikal  Rift  zone  and  the  Urals;  Figure  3). 

As  expected,  the  scattering  amplitudes,  A,  show  strong  variations  which  remarkably  well  correlate  with  tectonic 
areas  (Figure  3a).  High  A  values  are  found  in  the  Uralian  and  in  the  Sayan-Baikal  (labeled  BR  in  Figure  3b)  fold 
belts  which  could  potentially  be  related  to  stronger  surface  topography  and  complex  crustal  structures  due  to  folding. 
High  scattering  amplitudes  are  also  present  in  the  western  part  of  the  Siberian  Craton.  By  contrast,  the  eastern  part  of 
the  Siberian  Craton  (the  Minusisk-Aldan  High;  Figures  3a  and  b)  and  most  of  the  West  Siberian  Basin  show  low  A. 

The  geometrical  spreading  exponent,  y,  also  correlates  with  tectonics  (Figure  3b).  Note  that  the  level  of 
Yd  =  0.008  s"1  (purple  color  in  Figure  3b,  marked)  separates  the  stable  cratonic  areas  (with  y<  yD)  from  the 
tectonically  active  Baikal  Rift,  for  which  y  is  distinctly  higher  than  yD.  As  Morozov  et  al.  (in  press)  showed  by 
numerical  modeling,  such  a  difference  could  be  caused  by  generally  higher  level  of  contrasts  within  the  tectonically 
active  crust,  potentially  with  presence  of  relatively  thin  attenuative  zones.  As  these  authors  also  suggested,  within 
stable  parts  of  the  study  area,  y  varies  only  moderately,  from  about  0.004  to  0.007  s"1  (Figure  3b). 

The  effective  surface-consistent  coda  Qe  l  values,  defined  as  Qe~l  =  tdn ,  also  show  a  distinct  correlation  with 
tectonic  structures  and  reasonable  stability  within  the  study  area  (Figure  3c).  Coda  Qe  is  low  in  the  Baikal  Rift  zone 
(Qe~  100)  and  near  the  Urals  (Qe~  130)  but  higher  within  the  Siberian  Craton  (Qe~  300).  Within  the  eastern  part  of 
the  West  Siberian  Basin,  the  attenuation  is  the  lowest  ( Qe~  1000).  However,  although  indicating  a  reasonable 
correlation  with  the  geological  structures,  the  values  of  Qe~l  appear  somewhat  over-estimated  (could  be  to  a  factor  of 
about  2)  compared  to  detailed  measurements  performed  on  two  individual  PNEs  (Morozov  et  al.,  in  press).  Potential 
reasons  for  such  an  overestimation  will  be  considered  in  future  research. 
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Figure  3.  PNE  Lg  coda  Q  inversion  results:  a)  coda  scattering  amplitude  A;  b)  uncompensated  geometrical 
spreading  y,  2)  effective  coda  attenuation  Qj'  =  rdn:  In  plot  b),  note  the  level  of  yD  =  0.008  s'1  which 
was  suggested  to  discriminate  between  stable  and  active  tectonic  areas  (Morozov,  in  press).  Major 
tectonic  regions  are  indicated  in  plot  b):  BR  -  Baikal  Rift;  SC  -  Siberian  Craton;  MAH  -  Minusinsk- 
Aldan  High;  WSB  -  West  Siberian  Basin  PNE  profiles  are  labeled  in  plot  a). 


CONCLUSIONS  AND  RECOMMENDATIONS 

Morozov  et  al.  (in  press)  demonstrated  that  short-period  coda  attenuation  of  regional  arrivals  can  be  interpreted  as  a 
combination  of  uncompensated  geometrical  spreading  (y)  and  effective  attenuation  (Qe).  Taken  together,  these 
parameters  represent  the  generalized  attenuation  coefficient:  cdj)  =  y+  7iflQe,  which  in  many  cases  can  be  measured 
directly  from  the  data.  Parameter  ^removes  the  common  ambiguity  of  the  Q(f)  and  is  most  useful  for  correlating 
different  study  areas  and  tectonic  types.  Modeling  of  P-wave  coda  at  regional  distances  and  analysis  of  two  Peaceful 
Nuclear  Explosions  in  Russia  in  terms  of  this  concept  showed  a  remarkable  agreement,  especially  in  y  values. 


167 


2008  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


In  this  report,  we  applied  the  same  parameterization  to  the  inversion  of  spatially-variable  surface-consistent  Lg  coda 
parameters  from  five  PNE  profiles  in  Siberia.  The  resulting  map  of  the  geometrical  parameter  (;/)  is  in  excellent 
agreement  with  both  modeling  and  inversion  in  Part  I,  and  also  with  a  summary  of  over  forty  worldwide  studies 
(Morozov,  in  press).  In  addition,  scattering  amplitudes  were  found  to  be  highly  variable  and  strongly  correlate  with 
geological  structures.  Coda  Qe~l  values  also  showed  pronounced  correlation  with  geological  structures,  although 
these  values  also  appeared  somewhat  overestimated  compared  to  the  detailed  point  studies  in  Morozov  et  al.  (in 
press). 
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